Macroscopic quantum tunneling of two-component Bose-Einstein condensates 
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We show theoretically the existence of a metastable state and the possibility of decay to the 
ground state through macroscopic quantum tunneling in two-component Bose-Einstein condensates 
with repulsive interactions. Numerical analysis of the coupled Gross-Pitaevskii equations clarifies 
the metastable states whose configuration preserves or breaks the symmetry of the trapping poten- 
tial, depending on the interspecies interaction and the particle number. We calculate the tunneling 
decay rate of the metastable state by using the collective coordinate method under the WKB ap- 
proximation. Then the height of the energy barrier is estimated by the saddle point solution. It 
is found that macroscopic quantum tunneling is observable in a wide range of particle numbers. 
Macroscopic quantum coherence between two distinct states is discussed; this might give an addi- 
tional coherent property of two-component Rose condensed systems. Thermal effects on the decay 
rate are estimated. 

PACS numbers: 03.75.Fi, 05.30.Jp, 32.80.Pj 
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I. INTRODUCTION 

Multicomponent Bose-Einstein condensates (BECs) of 
alkali-metal atomic gases are expected to exhibit macro- 
scopic quantum phenomena that have not been found in 
a single condensate. Multicomponent atomic gases can 
be obtained experimentally by trapping different atomic 
species or the same atoms with different hyperfine spin 
states. The experimental realization of multicomponent 
BECs 1^, ^, ^ further stimulated many researchers to 
study the physics of this interesting system. 

Macroscopic quantum tunneling (MQT) is an inter- 
esting subject in many fields of physics. In this paper 
we study MQT of metastable two-component BECs in a 
trapping potential. Thus we need to know detailed in- 
formation about the stationary state of this system. The 
structure of the ground state has been studied by solving 
two coupled Gross-Pitaevskii equations (GPEs) analyti- 
cally or numerically g 1 1, 0, |[ |, , |l|, g |l|, ^ . 
The stationary solution of the GPEs gives the density 
profile of the condensate characterized by the parame- 
ters of the system — trapping frequencies, the number of 
atoms of each component, and three s-wave scattering 
lengths, ai, a2, and ai2, which represent the interactions 
between like and unlike components. 

The interspecies interaction characterized by ai2 plays 
an important role in determining the structure of the 
ground state. When the inequality ai2 > ^JoiOa is satis- 
fied, a mixture of two-component BECs without a trap- 
ping potential tends to separate spatially The 
trapped BECs have two different configurations of the 
condensates when ai2 is large 0. One configuration 
preserves the spatial symmetry of the trapping potential 
by forming a core-shell structure. The other breaks the 
spatial symmetry by displacing the center of each con- 
densate from that of the trapping potential. 

Ho and Shenoy first constructed a simple algorithm to 



determine the density profile within the Thomas-Fermi 
approximation (TEA) g. However, the TFA is not 
enough to describe the density profile of phase separa- 
tion because the penetration at the boundary of each 
component is not considered. Without the TFA Pu and 
Bigelow investigated numerically the ground state of a 
Rb-Na BECs by assuming spherical symmetry [|| . When 
ai2 is large, they found a ground state that forms a core 
of Rb at the center of the trap and a shell of Na around 
Rb, and a metastable state that has a Rb shell and Na 
core. However, they noted the existence of an unstable 
mode which forms the core-shell structure. After that, 
further investigation of two- or three-dimensional GPEs 
showed a spherical symmetry-breaking solution for the 
true ground state 0, |l^, 

Ohberg showed that whether the ground state takes 
a symmetry-breaking state (SBS) or a symmetry- 
preserving state (SPS) depends not only on the inter- 
species interaction but also on the particle number, the 
intraspecies interaction, and the shape of the trapping 
potential ||l^. However, the detaila of the metastable 
state have not been studied. Thus, we investigate the de- 
pendence of the ground state and the metastable state of 
two-component BECs in a cigar-shaped potential, which 
can be considered as a quasi-one-dimensional system for 
simplicity. We also make a linear stability analysis of 
the stationary solutions of the GPEs and reveal their 
metastability. 

A metastable BEG can also be found in a single con- 
densate with negative s-wave scattering length [|l^ . The 
negative scattering length represents an attractive atom- 
atom interaction, which causes the condensate to collapse 
upon itself to a denser phase. The balance between the 
attractive interaction energy and the zero-point kinetic 
energy of the trapping potential realizes the metastable 
condensate. MQT of a condensate with attractive inter- 
action has been predicted |jl8[. 
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For two-component BECs with repulsive interactions 
the metastabiUty mainly comes from the competition be- 
tween intra- and interspecies interactions. We study the 
transition between the SBS and SPS by MQT. 

In Sec. II, we obtain the stationary solution of the 
GPEs numerically. The phase diagram of the ground 
state has a rich structure including metastable states. 
The stability of these solutions is checked by following 
Ref. [0 which considers the stability by taking account 
of the linear fluctuation around the stationary solution. 

In 



with the chemical potential fj,i . Here the two-body inter- 
actions Uii and Ui2 are written as 



Uii = gu / dydz\^ij_{y,z)\ 



9ii 
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(3a) 



Ui2 = 5i2 / dydz\il:i^{y,z)\'^ \'ip2±{y,z)\'^ 

912 



(3b) 



Sec. Q we introduce the collective coordinate where bi± = y/h/miUj^±, g^i = Airh^ai/rai and 512 



method to evaluate the decay rate of a metastable state 
through MQT by an imaginary-time path integral (in- 
stanton) technique. The collective coordinate enables 
us to derive an effective one-dimensional Lagrangian de- 
scribing the two-component BECs and obtain the decay 
rate. We estimate the decay rate at finite temperatures; 
the results show the probability of observation of MQT. 
We also discuss macroscopic quantum coherence (MQC), 
which is the oscillation between the SPS and the SBS. 



Section [V is devoted to conclusions and discussion. 



II. FORMULATION AND STATIONARY 
SOLUTION 



2Tih ai2/mi2 with reduced mass TO12. The corresponding 
two coupled time-dependent GPEs are given by 



(4) 



where 
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-ff2[V'l:'02] = - 



2m2 dx"^ 2 



+ -'m2Ujlx^ - fl2 



+ C/22|l/'2p + C/l2|V'l 



(5a) 



(5b) 



We consider two-component BECs in the external 
trapping potentials 



and each wave function is normalized by the number of 
particles Ni as / dx\'ipi{x)\'^ = Ni. 



1 



1 



Kxt(r) = :^m,iufx' + -m^LufJy' + z'), z = 1, 2, (1) 



Numerical solution 



where rrii is the atomic mass, and uji and uji± are the 
longitudinal and transverse trapping frequencies. For 
u!i± ^ uji the trapping potential is cigar-shaped. If 
the two-body interaction energy is smaller than huj±, 
it does not affect the transverse component '>pi± of the 
wave functions, which allows us to analyze the prob- 
lem in one-dimensional space. Although it has been 
predicted that the two-body interaction is changed by 
the effect of tight confinement of the trapping potential 
0, we will use the following treatment to derive 
the onc-dimensional GPEs Using the ground state 
wave function in the harmonic potential for ^i^{y,z), 
we assume the macroscopic wave function as 4'i(r,t) — 
ipi(x, t)ipi±(y, z) {i = 1, 2). These wave functions are sub- 
stituted into the three-dimensional Gross-Pitaevskii en- 
ergy functional, which is integrated over y and z. Thus 
we obtain the one-dimensional Gross-Pitaevskii energy 
functional, 



'^[V'1,'02] = / dx 



^ \ 2to,- 



dx 



, 1 2 2| ; |2 



-M^IV'^n +^C/ll|V'l|' + ^f/22|^2|* 



+ Ui2\^Pl\^\^2 



(2) 



The stationary solutions of Eq. (|j) correspond to the 
critical points of the energy functional H. There are sev- 
eral ways to find these critical points numerically. Our 
method is described in the following. The stationary so- 
lution ipio satisfies the relation 



^^-i[V'lO,'02o]'0jO = 



(6) 



from Eq. (|j). The solution is taken to be real by 
making the phase zero. Using the trial function ipl", 
Tpio is given by ipl" and its deviation Aipi, i.e., jpio = 
^tn — Alpi. Substituting this relation in Eq. (^, we 
obtain the linearized equation for Aipi: 



ffi[<',V'r]+2{/ii(v^rTAV'i 

-H2C/i2VfV2"Al/'2 =ai, 
i?2[^r,V'2"]+2C/22(V'2")'AV'2 
+2C/i2VfV2"A^l =(72, 



(7a) 
(7b) 



where cr, = Hi[tpl" , ■ip2"]i^i" ■ The linear correction Atpi 
can easily be calculated and the modified trial function 
is defined by ipi — tpl" — Aipi. We repeat the above 
calculation until the solution converges by conserving the 
norm of each component. 

Assuming the condensates of two hyperfine spin states 
of ^^Rb, we use the values of the scattering lengths 
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ai = 5.36 nm and 02 = 5.66 nm, which have the ra- 
tio a-ija-y — 1.06 We choose the atomic mass 
TOi = TO2 = TORb = 1-45 X 10"^^ kg, the trapping fre- 
quency wi = ti;2 = tf = 90 X 27r Hz, and the aspect ratio 
lo^/lo = 30. It is convenient to introduce scales char- 
acterizing the trapping potential: (a) the length scale 
b = yj h / mi2Uj , (b) the time scale uj~^, and (c) the en- 
ergy scale huj. The dimensionless parameters normalized 
by these scales are expressed by putting a tilde upon the 
symbols. Then the dimensionless intraspecies interac- 
tions become Uu = 0.2010 and U22 = 0.2123 from Eqs. 
(^. Setting the particle numbers to A^i = A^2 = for 
simplicity, our formulation has two free parameters, N 
and Ui2. The parameter U12 might be controlled exper- 
imentally by the choice of some combination of atoms, 
or by changing the scattering length via the Feshbach 
resonance ||2^. 




FIG. 1: Numerical solutions of the two coupled stationary 
GPEs. The solid lines and the dotted lines show the conden- 
sates 1 and 2, respectively. The wave functions ipi divided by 
y/N and the coordinate x are normalized by the length scale 
h = --Jhlnwiuj. (c) shows the symmetry-preserving state with 
energy 87.4261 in units of TiujN . (d) shows the symmetry- 
breaking state with energy 87.5514. 

Typical stationary solutions of Eq. (Q) are shown in 
Fig. l[ When the interspecies interaction is weak, two 
condensates overlap each other. For U12 < Uu < U22 
two overlapping condensates have peaks of the density 
at the center of the trapping potential, as shown in Fig. 

|l|(a). For till < U12 < U22 and \/UiiU22 < U12, the 
density peak of condensate 2 is not at the center and 
two peaks appear symmetrically about the origin x — Q, 
as shown in Fig. |l|(b). Note that the width of conden- 
sate 2 is larger than that of 1, because JJu < U22- These 
structures can be predicted easily within the TFA |Q, ^ . 

For C/12 > \/lJ\\lJ22 the two condensates separate from 
each other with very narrow overlapping regions. In Fig. 
|l|(c) the condensate 1 occupies the central region, pushing 
aside the condensate 2 symmetrically; this configuration 
preserves the spatial symmetry of the trapping poten- 



tial. On the other hand, as shown in Fig. |^(d), there 
exists another configuration with the boundary between 
the two condensates at the center of the trapping poten- 
tial and its spatial symmetry is broken. Now we call the 
stationary state in Fig. 1(c) the symmetry-preserving 
state and that in Fig. ^d) the symmetry-breaking state. 
The total energy of solution (c) is lower than that of (d) 
as described in the figure caption, so that the solution 
(d) represents the metastable state. 
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FIG. 2: N-U\2 phase diagram of the ground state. The 
gray region represents the overlapping configuration and the 
other the separated configuration. In region B the SBS is the 
ground state, while in region P the SPS is the ground state. 
In the region with slanted fines, there exists the metastable 
SBS (SPS) denoted by the lower-case letter b{jp). 

In Fig. H we show the N-U12 phase diagram of the 
ground state. The gray region represents the overlapping 
configurations. Fig. ^(a) and|l|(b), and the other the sep- 
arated configurations. Fig. |l|(c) a nd [l|(d) ; the two regions 

are divided by the hne U12 — VU11U22 = 0.2066 which 
was predicted by the TFA §, |, § 0. The gray 
region is further divided into two regions [Fig. ^(a) and 
|l|(b)] by the line C/12 = 0.2010 ||l|. More precisely, these 
boundaries are bent for small N because of the failure of 
the TFA. The region of the separated configurations has 
the following structure. The bold line shows the bound- 
ary where the energy of the SBS is equal to that of the 
SPS. In region B the SBS is the ground state, while in 
region P the SPS is the ground state. 

The position of the SBS and the SPS in the phase di- 
agram Fig. 1^ can be understood as follows. We first 
consider the transition on increasing the particle number 
N with a fixed value of U12 ■ Note that the SBS has one 
domain wall and the SPS two domain walls. When N is 
small, the SBS is realized because the multiple domain 
walls increase the domain wall energy, which is estimated 
by the energy / rfa:E,(?iV2mi)|VV^P + C/i2|^iP|V'2p] in 
Eq. (H). The increase in N makes the intraspecies in- 
teraction energy important, thus tending to extend each 
domain. This overcomes the energy of formation of do- 
main walls, so that the SPS becomes more stable than 
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FIG. 3: Schematic illustration of the total energy in the region 
of the separated configurations in Fig. ^. When the energy of 
the SBS is equal to that of the SPS, the energy configuration 
becomes a triple well. 



the SBS. When U12 increases, the domain wall energy 
becomes larger and thus the region B is extended. 

The bold line suggests the existence of metastable 
states as shown in Fig. ^. In Fig. ^, the regions Bp 
and bP have metastable states; the SBS is the ground 
state and the SPS is the metastable state in the region 
Bp, and vice versa in the region bP. The details of the 
metastable state and how to decide the boundaries be- 
tween B and Bp, bP and P are described in the next 
subsection. 



B. Stability of the solutions 

We linearize the energy functional Ti of Eq. (^) by 
substituting 



-01 = + Hi- 



(8) 



Here ipio is the stationary solution obtained by solving 
the GPEs, and the fluctuation JV'i is complex. The sta- 
tionary solutions represent the local minima or the saddle 
points of the energy functional. Then the energy can be 



expanded around the stationary solution: 

n[^i,^2] ^ no + 6n[s^bi,s^P2l (9) 
1 r 

Here -q (f^i, = {5il)i,5il)2,Hl,H2) and 

W = {Wij) is the Hessian operator, corresponding to 
the second-order derivative of the energy functional at 
the stationary solution: 



Wii 
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W22 = W44 = 



2mi dx^ 



1 22 



Ml 



-f2i7iiV^?0 + C^12V'20> 



1 



M2 



1^12 
1^13 



+ '2U224>lo + Ui2i)l^, 
M^23 = VF34 = VF41 = C/l2V'10V'20, 



10' 



W2A ^U22i>: 



20- 



When all eigenvalues of W are positive, the stationary so- 
lution is stable, while the appearance of a negative eigen- 
value makes it unstable. 

This eigenvalue problem is simplified by the unitary 
transformation [|16| 



u'fwu 



Li 
L2 



where 



L2 



Wn - C/iiV'' 
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W22 - C/22V'20 



T^ii + t/iiV'?o 2W^i2 

214^21 M^22 + U22^4o 



(11) 

(12) 
(13) 



and 



U 



a b 
a b 
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I) 



Law et al. used the lowest eigenvalue of L2 as the stability 
criterion of the system of two-component BECs |l^ . The 
lowest eigenvalue of Li is zero and the eigenfunction is 
given by the stationary solution ipiQ. 

Figure ^ shows several lower eigenvalues of W as func- 
tions of TV for the SPS (a) and SBS (b) with U12 = 0.2438 
as used in Fig. |l|(c) and 0(d). The critical particle num- 
ber Nc defined by the zero eigenvalue of L2 gives the cri- 
terion for the stability of the stationary solution. From 
Fig. 1(a), the SPS is stable for N > N^. Figure |(b) 
shows that there always exists one negative eigenvalue 
whose eigenfunction changes N. Hence, as long as N 
is fixed, the SBS is stable for N < Nc. Obtaining Nc 
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as a function of U12 allows us to decide the boundaries 
between the regions B and Bp, P and bP in Fig. |. 

Finally, let us note the fluctuation changing the parti- 
cle number. By using the eigenfunction u — (7/1,1/2) of 
L2 this fluctuation is evaluated as 



■>Jjia{x)ui{x)dx. 



(14) 



For the mode in Fig. ^(b) whose eigenvalue is always 
negative we obtain SNi 7^ 0. The other mode of L2 in Fig. 
^(b) conserves the particle number. The fluctuation that 
changes the particle number leads to the ground state of 
the SBS with unbalanced particle number A^i 7^ A''2 [M. 



III. POSSIBILITY OF MACROSCOPIC 
QUANTUM TUNNELING 

As described in Sec. II, the SBS is the ground state 
and the SPS is the metastable state in the region Bp in 
Fig. 1^, and vice versa in the region bP. In this section, 
we study the MQT of the metastable state in Bp and bP. 
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FIG. 4: Several lower eigenvalues of W for the SPS (a) and the 
SBS (b). The bold lines show the eigenvalues of L2 and the 
dashed lines show those of Li. The critical particle number 
is represented by Nc. 



A. Collective coordinate approach 

It is difficult to consider MQT by full quantum field 
theory. In the case of a single condensate, the variational 
method is often used to estimate the condensate wave 
function. This method was applied to the evaluation of 
the decay rate via MQT of a metastable condensate with 
attractive interaction [ p^ . However, there is an impor- 
tant difference in the description of the decay of a sin- 
gle condensate and the transition between the SBS and 
SPS in two-component condensates. In the former case, 
there is an obvious collective coordinate, i.e., the spa- 
tial size of the condensate wave function, which allows 
us to approximate the wave function under the Gaussian 
ansatz |l^. In contrast, in the latter case, it is difficult 
to find suitable collective coordinates that can describe 
the continuous deformation from a metastable state to 
the ground state. Thus we introduce an alternative vari- 
ational approach for this system in order to calculate the 
MQT rate. 

The action S for the Gross-Pitaevskii model is given 
by 5* = / Cdt with the Lagrangian 



i=l,2 



dx 



(15) 



where Ti. is the Hamiltonian of Eq. (^ . The macroscopic 
wave functions are written as tpi = <&i(a;, t)e'^*'^'*'', where 
<I>f = Pi and 9i are the number density and the phase for 



each component, respectively, 
in Eq. (p|) yields 



Substituting these forms 



j dx 






1=1.2 K 



2m, 



9^ 

dx 



where V{^i) is written as 



(16) 



+ ^$4 + ^$4^t/i2$?<i>i (17) 

The amplitude $i can be expanded around the sta- 
tionary solution = \/Pi ~ IV'iol in Sec. II by using an 
orthogonal complete set Uin{x), 

^,{x,t) ^ ^l(x)+^Qn{t)u,n{.x) (n=l,2,...), (18) 

n 

with the normalization 

^ ^ / '^in^iradx — (^nm- (-i^) 
i=l,2 •' 

Here Qn{t) stands for the dimensionless arbitrary func- 
tion and represents the small displacement of the density 
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profile from the stationary solution: 

E^nW= E / dx[Mx,t)-^U^)f (20) 

(n= 1,2,...). 



i=l,2 ' 



Substituting Eq. ( [l8|) to Eq. ([l6|), we can obtain 
the effective action written by the functions Qn{t)- If 
we assume that the phase has the form 9i{x,t) 
12n^ri{t)vin{x)/N with somc complete set Vin{x), the 
first term of Eq. ( p^ can be written as 

Y^2n^^Qn{t)J2 J dxv,M<^t{x)u„,{x), (21) 



by using pi ~ 2 (5„$^iti„. Then we define the col- 
lective coordinate Qn = bQn/VN and the collective mo- 
mentum Pm, = 2UPm/b with the length scale 6 of the 
trapping potential. By choosing the complete set Vin(x) 
as 



> ^ y dxVijn(X) ^ Uin{ 



x) = <5„„, (22) 



the first term of Eq. ( |l6| ) becomes PnQn and the 
second term 



Prn(^)Pn(t) 
9 A4 ' 

^iy± rn.7r 



(23) 



where the effective mass matrix Mmn is given by 
62 



1—1,2 



g I dVirn dVifi 



dx. (24) 



The potential V{^i) is a function of the collective coor- 
dinate Qn- Thus we can obtain the effective action 



dt 



E Pr^Qn - Wcff(P, Q) 

n 

P = (Pi,P2,...), Q = (Qi,Q2,...), 

with the effective Hamiltonian 

P P 



(25) 



Hefl(P,Q)=E 



2Af„ 



y(Q). 



(26) 



By substituting Eq. (|T8|) into Eq. (|T^ the effective 
potential ^(Q) can be expanded as follows: 

nQ) = / m)dx+iE^«E"™^^^"^"+---- (27) 



Here the linear term in Qn vanishes because $| is the 
stationary solution. The quadratic term of Q„ can be 
written by the Hessian operator TJ^j, which is equal to 
2L2 given by Eq. (p^). We take the orthogonal set Uin 



as the eigenfunction of L^. Thus the second term of Eq. 
( ^ ) is diagonalized and T^(Q) is written as 



nQ) = / m)rf^ + iE^^"^ 



(28) 



where A^j in = 1,2, ...) are the eigenvalues of Hij. The 
constant / V{^l)dx will be chosen to be zero in the fol- 
lowing section. 



B. Calculation of the decay rate 

We now calculate the MQT rate F of the metastable 
state in Fig. |. The energy of the metastable state must 
have an (exponentially small) imaginary part when the 
tunneling is taken into account. Then the decay rate of 
the metastable state is given by [p4l 



F . -Imi^,. 



(29) 



The energy Eg is evaluated by the partition function 



as 



lim 

(3^00 



/3 ' 



(30) 



(31) 



where f3 — 1/kBT. Using the action Eq. ( |25| ) with the 
imaginary time t —^ —it (Euclidean action ^b), the par- 
tition function is written as 



Z{(i) 



j DQ{t) exp 



SeIQ) 



(32) 



Within the WKB approximation this path integral is 
evaluated by the saddle-point method j2^. More pre- 
cisely, the dominant contributions to the path integral 
are from paths that minimize the Euclidean action Se- 
Such paths are the solution of the Euler-Lagrange equa- 
tion 6Se/SQ = 0, the classical equation of motion for 
the valuables Q(r) in the inverted potential — V^(Q). We 
choose the boundary condition that Q(t) approaches the 
metastable minimum at r = ±00. By solving this equa- 
tion of motion, we obtain the solution Qb called the 
"bounce solution." Then the decay rate has the form 



F ~ ^ exp 



Sb 



(33) 



where Sb = Se{Qb) is the Euclidean action evaluated 
at the bounce solution Qs and A the quadratic quantum 
fluctuation around the bounce solution. The following 
describes how to approximate the bounce solution. 

We are interested in the regions Bp and bP near the 
dashed lines in Fig. |[ These regions have metastable 
states as described in Sec. H. On the dashed lines, one 
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of the eigenvalues of the Hessian operator Hij vanishes. 
Thus, in the region close to the dashed lines, the potential 
barrier of the metastable state is very small along the di- 
rection of the eigenfunction with the zero eigenvalue. We 
may assume that the direction of the initial (infinitesi- 
mal) velocity of the bounce solution is given by the eigen- 
function Ui subject to the following conditions. First, Ui 
has the eigenvalue which is small in this region and 
becomes zero on the dashed line. Secondly, this eigen- 
function conserves the particle number, i.e., SNi = in 
the sense of Eq. (p^. Thus the trajectory of the bounce 
solution is mainly described by the collective coordinate 
Qi{t), which is the coordinate along the direction of Ui, 
and the other coordinates Q2{t),Q3{t), . . . give higher 
order corrections of the solution. These assumptions al- 
low us to solve the bounce solution approximately; the 
trajectory of the bounce solution is straight in the collec- 
tive coordinate space Q = (Qi, (52, Qs- ■ • ) [psl. Thus 



the infinite-dimensional system of Eq. (25) is reduced to 
a one-dimensional quantum mechanical system with the 
collective coordinate Qi subject to the action 



dt 



M fdQi 



V{Qi 



(34) 



Here we have defined the mass M by the (1,1) component 
of Eq. (^J); the other components represent the mass 
relevant to Q2,Q3 . . . , so that they are negligible. 
The mass M includes unknown functions vn, although 
they satisfy Eq. (p2[). We will assume vn ~ 0(1), thus 
obtaining M = Mn - {AmmN'^ /b'^) x (b'^/N) - niRbN. 

The next problem is to decide the form of the effec- 
tive potential V{Qi). It should be noticed the poten- 
tial is expressed as ^(Qi) ~ ^NXKQl/b'^toi small Qi 
around the metastable state, from Eq. (p^). As Qi in- 
creases, it is not clear how to extrapolate the potential. 
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FIG. 5: The quadratic-plus-cubic potential given by Eq. (p5[). 
The potential has a metastable minimum at Qi = 0, barrier 
height Uo = at Qi = Qf, and width Rq. The second 
derivative of V{Qi) at Qi = is given by XiN. The coordi- 
nate Qi corresponds to the sphaleron with energy AV^. 
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FIG. 6: Saddle-point solution (sphaleron) of the two coupled 
stationary GPEs for U12 = 0.2438 and TV = 2000. The unit 
of length is the same as in Fig.l. The inset shows the detail 
near x — 8. 



However, V{Qi) should reflect the structure of the origi- 
nal potential V{^i) of Eq. ([T7|), which has a metastable 
state in addition to the ground state as discussed in Sec. 
11. Hence we require, first, that V{Qi) has a metastable 
minimum at Qi = 0. Secondly, as Qi increase, V{Qi) in- 
creases once and decreases via a potential barrier Ay as 
shown in Fig. ||. The potential is expected to be written 
as a power series of Qi- Since the calculation of the de- 
cay rate does not need information on the ground state, 
we neglect the nth order terms (n > 4) and approximate 
the potential as 



ViQi) 



-NXI 



1 
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aQl (a>0). 



(35) 



with an unknown parameter a. Then the height of the 
potential barrier is given by 



AV 



(36) 



The value of AV^ cannot be determined within the col- 
lective coordinate method. The barrier AV^ may be in- 
terpreted as follows. In general, when we consider quan- 
tum tunneling, it is natural to assume that the bounce 
trajectory will go through the sphaleron, which is the 
unstable stationary solution of the equation of motion, 
corresponding to the saddle point of the potential [ p6[ . 
Then AV represents the energy of the sphaleron. In our 
case, the equation of motion is the GPE of the poten- 
tial Eq. (|l7|); we have found the sphaleron by numerical 
simulations (see Fig. ||) and obtained the value of AV. 
The collective coordinate Qi of the sphaleron is simply 
written as 



n = 2^ 



/ 3AV^ 
2XlN' 



(37) 



Thus our collective coordinate Qi effectively describe 
MQT: the points Qi = and Qi = Qf correspond to 
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the metastable state and the sphaleron, respectively, and 
the tunneHng is represented by the bounce solution, going 
through the sphaleron. We will give the explicit bounce 
solution written via the collective coordinate in Eq. (44). 

To calculate the decay rate of Eq. (^3|), it is conve- 
nient to introduce new scales characterizing the quantum 
tunneling instead of the scales of the trapping potential: 
according to Fig. [s] we define the length scale 



I 3AV 



(38) 



the energy scale 

Uo = AV, (39) 
and a time scale representing the "tunneling time" 



To 




mi2iV 



(40) 



with ujQ = y/27/2. In these units the action Eq. ( |34| ) can 
be written by the dimensionless length q = Qi/Rq and 
the time s — I/tq as 



S 

h 



ds 



Viq) = -uj',q\l-q). 
Here h is the effective Planck constant defined as 

n 



h 



ToUo 



1 

UJQ 



mi2N ^yhLLl\l 



M 



AV 



(41) 
(42) 

(43) 



whose value must be smaller than unity for use of the 
WKB approximation, although it includes the macro- 
scopic valuable N. From Eq. ( |4l| ) the bounce solution 
can easily be obtained by solving the equation of motion 
Cpq/ds^ = dV{q)/dq with the boundary condition q~Q 
at s = ±00: 



qB{s) = sech^ ' 
and the decay rate can be written as 



A 



To 



■ exp 




(44) 



(45) 



with the prefactor A — Ay^ajQ/nh and the bounce action 
Sb = 8aJo/15- 

We may observe MQT experimentally if the decay rate 
F is of the order of the lifetime of the BEC. Let us search 
the region near the dashed lines in Fig.|| satisfying this 
condition. We obtained the potential barrier AV from 
the sphaleron energy and the eigenvalue from the Hes- 
sian operator Hij . Recalling that N is equal to the criti- 
cal particle number Nc on the dashed line (Sec. II B), it is 



convenient to introduce a small parameter i5 = |1— A''/7Vc|. 
Figure shows the 5 dependence of Af and AV for the 
metastable SPS and SBS near the dashed lines in Fig. ||. 
Since Af and AV vanish on the dashed line with S — 0, 
the scaling laws for the particle number are expected to 
be like those of a single condensate E^: 



A? 
AV 



(46) 
(47) 



The exponents /3,7 and the coefficients C,V are deter- 
mined by fitting the scaling laws to the numerical re- 
sults. Thus we obtain the exponents (3 = 1.0±0.002, 7 — 
2.06±0.03 for the metastable SBS, and (5 = 0.788±0.006, 
7 = 1.673±0.007 for the metastable SPS. As shown in Ta- 
ble. |, they are approximately independent of the value 
of U12 within our analysis, while the coefficients C and 
V depend on C/12. When we calculate i?o and Uq of Eqs. 
( p8| ) and ( p9| ) by using these exponents and coefficients, 
we find that the metastable SPS has larger values of i?o 
and Uq than the SBS. Thus MQT cannot be expected for 
the metastable SPS compared with the metastable SBS. 

Substituting Eq. (|46|) and Eq. (^) to tq and h of Eqs. 
( ^ ) and ( p3| ) , we can obtain the scaling laws of the decay 
rate F of Eq. ^ by 



Sb 



8c^o 



M V 



15 V rni2N ^ 



and 



A 

To 



4wn 



M 



TT V TOl2iV 



1/4 



|^^p2^ 1/4^/3/4+7/2^ 



(48) 



(49) 



The dominant contribution to F is the exponential 
factor. To obtain the observable decay rate, we re- 
quire h ~ 10^^, although h should be small under 





0.1 

s 



0.1 



FIG. 7: 5 dependence of the eigenvalue Ai and the potential 
barrier AV . The metastable SPS is shown in the left column 
and the metastable SBS in the right column. The open circles 
represent the numerical results of Ai and AT/ from the GPEs. 
The solid, dashed-dotted, and dotted lines show the scaling 
laws for i/12 = 0.2250, 0.2483, and 0.2813, respectively. 
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TABLE I: The values C, V, 13, andj of Eqs. @ and @ 
metastable SPS 



Ul2 


0.2250 


0.2483 


0.2813 


c 


0.03041 


0.021508 


0.015738 


V 


57.155 


139.87 


279.79 


p 


0.78718 


0.79373 


0.79177 


1 


1.6684 


1.6801 


1.671 


metastable SBS 




0.2250 


0.2483 


0.2813 


C 


1.1445 


2.0161 


3.5006 


V 


17.454 


20.462 


26.086 


P 


1.0034 


1.0022 


1.0055 


7 


2.068 


2.0371 


2.0562 




0.01 . \ 



(a) 



C. Macroscopic quantum coherence 

An interesting phenomenon may appear on the bold 
line in Fig. || where the energy of the SPS is equal to 
that of the SBS. It is macroscopic quantum coherence be- 
tween the SBS and the SPS, i.e., the oscillation of a wave 
packet between their potential wells. Then the effective 
potential V{Q) has a triple-well geometry, as illustrated 
in Fig. |. 

The period of that oscillation can be estimated by the 
splitting A of the ground state energy due to the tun- 
neling. The splitting is written as A ~ Ae~^^^ , where h 
is the effective Planck constant, / the instanton action, 
and A a prefactor of the order of unity. Here we only 
estimate the order of the oscillation period by using the 

ofEq. (11). 



effective Planck constant h — fi/ {R^f MUq) 

The barrier height Uq is given by the energy of the 
saddle-point solution, and the barrier width R is given 
by the "distance" between two stable solutions of the 
SBS and the SPS. Then the distance is estimated to be 
R = {Ri + i?2)/2 as shown in Fig. ||, where Ri and R2 
are given by Eq. (38). By tuning the parameter N 200 
and U12 ~ 0.2116, the period of the oscillation becomes 
of the order of 1 sec. The increase of the particle number 
raises the energy barrier between the SPS and the SBS 
so that MQC cannot be observed. 



0.001 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 

S =|l-N/N 



D. Finite-temperature effect 



FIG. 8: The solid lines show the effective Planck constant of 
the metastable SBS and the dashed lines show that of the 
metastable SPS. The parameter U12 is set as 0.2250 for (a) 
and (d), 0.2438 for (b) and (e), and 0.2813 for (c) and (f). 
The critical particle number Nc is (a) 481, (b) 770, (c) 1197, 
(d) 1121, (e) 2520, and (f) 5575. In the region h > 1, the 
WKB approximation breaks down. 



Let us consider finite-temperature effects, although the 
discussion of the last subsection was limited to zero tem- 
perature. Then the bounce solution Eq. (44) turns into 
the periodic solution, i.e., the classical solution in the po- 
tential -V{q) [Eq. (H)] with energy -E {0 < E < 1) 
[pTf. From Fig. Il^ the explicit solution is given by the 



the WKB approximation. Figure shows the effec- 
tive Planck constant for several values of U12 as a func- 
tion of 5 = |1 - N/Nc\. The SBS has a wider range 
with respect to 6 satisfying the above condition for h 
than the SPS. Although it is difficult to tunc the value 
of 6 experimentally, we may observe MQT of the SBS 
more easily than that of the SPS. We now estimate the 
range of S for U12 — 0.2483 where P becomes of the 
order of 10"^ sec^^; the lifetime of the BEC is typi- 
cally 100 sec 0. The metastable SPS {Nc = 770) yields 
Ss/h = 9713 X (51-283 and A/to = 44.66 x w x 5^ °^^ sec~\ 
so that 6 < 3.9 x lO^'^, a range too narrow to observe 
MQT. For the metastable SBS {Nc = 2520) we obtain 
Ss/h = 146.7 X 51-537 and A/tq = 53.15 x w x Ji-^es 
sec~i; then S < 0.205. However, values of (5 < 6.4 x 10^^ 
make h larger than unity, thus breaking the WKB ap- 
proximation. 



Energy 

















A V 





SBS u R SHS 

^ ^ SPS 



2 2 

FIG. 9: Schematic illustration of the triple-well potential. A 
shows the splitting of the ground state energy. See the text 
about Ri and R2. 
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elliptic function 

= 92 - (92 ~ 9i)sn2 



9o.s; m 



(50) 



with m = \/ (92 ~ Qi)/{<12 — 9o)i and the period T is given 
by the complete elliptic integral of the first kind K{m), 



K{m). 



(51) 



The period T is related to the temperature, T — hf3, 
where (3 is the inverse temperature normalized by Uq. 
This solution reduces, of course, to the previous solution 
Eq. (Q) for E — 0. The corresponding bounce action is 
evaluated as 



SbW = / ds 
Jo 



ds 



= W + hpE, 



(52) 
(53) 



where 
W 



15 



W0V92 



90 



X [2(9o + 9i + 92 - 9091 - 9o92 - qiq2)E{m) 
+ (91 - 9o)(29o - 91 - q2)K{m)] (54) 

with the complete elliptic integral of the second kind 
E{m). The decay rate by MQT takes the form 



The prefactor A{[3) was derived in Ref. | 



(55) 



as 



A(/3) 



£|(92 - 90)^/^(92 -9i)(l-™^) 



X [a{m)E{m) + b{m)K{m)]-'^'^ sinh (^^^^ j (56) 
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FIG. 10: Turning points qi and 92 in the potential —V{q) — 
-(l/2)a;g(?2(l - q), loo = ^/27/2. The "energy" E {0 < E < 
1) is determined as a function of inverse temperature fH by 
requiring that the motion between the turning points qi and 
q2 is periodic, with period f3h. 



with 



a{m) = 2(m^ - + 1), 
b{m) = {l-m^){m^ - 1). 



(57) 
(58) 



The thermal effect increases the MQT rate by a factor 
of only the order of unity from the MQT rate at zero 
temperature. 

For £; ^ 0, we have (1 - m^) sinh(wo/3V2) ^ 8, 
a(jn)E{m) + b{m)K{m) 2, and 90,91 ^ 0, 92 ^ 1, 
so that A{(3) A^Joj^j-Kh, which reproduces the zero- 
temperature decay rate T of Eq 
limit i? — > 1, where the period be 



Let us turn to the 
tiaves as 



2n 



(59) 



The leading term gives the crossover temperature f3^^ = 
hujo/2n. As the temperature is raised above /3~^, 
the system has no bounce solutions, and the decay is 
caused by thermal activation (the Arrhenius law): Tp ~ 
Wo exp(— /3AT^). The crossover temperature is of the or- 
der of O.lnK for the range of S discussed in the last sub- 
section. 

Equation ( |55| ) and the Arrhenius formula arc not avail- 
able in the narrow region near Pc- In this region the decay 
rate is given by ]29| 



r(/3)ro 



15/i7r2 



X exp 



V 2 

18/3c 




Pc) 



with the error function 
erf(a;) 



/27r J- 



dy exp 



' 2 



(60) 



(61) 



For small h <C 10^^, this formula matches smoothly onto 
Eq. (^) for P > l3c and the Arrhenius formula for /3 <_0c 
near /3c- However, we cannot apply the formula Eq. (|60|) 
to MQT since the value of h in our situation is of the 
order of 10~^. We leave the issue of the crossover region 
for future study. 



IV. CONCLUSIONS AND DISCUSSION 

The metastability and MQT of two-component BECs 
were studied theoretically. By analyzing two coupled 
GPEs numerically, we obtained two kinds of metastable 
state, the symmetry-breaking state (SBS) and the 
symmetry-preserving states (SPS), which depend on the 
particle numbers and the interspecies interaction. We in- 
troduced the collective coordinate method by improving 
the usual Gaussian variational approach, and calculated 
the MQT rate within the WKB approximation. The ef- 
fective potential V{Q) was determined by analysis of the 
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linear stability and using the saddle-point solution. Then 
the decay rate is found to obey a scaling law near the crit- 
ical region. MQT from the SBS to the SPS is expected to 
be observed in a wide range of the parameter 5. We also 
predicted MQC between the SBS and the SPS, although 
the range of 5 is rather narrow. 

Our analysis is restricted to the one-dimensional con- 
densate, but it can be applied to a system in a highly 
anisotropic trapping potential. The extension to the 
three-dimensional system is troublesome. However, the 
qualitative nature will be the same as in the one- 
dimensional case. This analysis can also be applied to 
the MQT between the domains of a spinor condensate 
where the external magnetic field can be used as an- 
other variable parameter. 

In Sec. II B we stated that negative eigenvalues always 
exist for the SBS, corresponding to the change of particle 
number. This instability may be caused by inelastic col- 
lisions of atoms in a real system. However, if we confine 



ourselves to the region near the critical particle number, 
MQT is expected to be the dominant mechanism of de- 
cay ||l^ . Thus the change of particle number is neglected 
in the analysis of the MQT. 

Finally, we comment on the validity of the quasi-one- 
dimensional approximation. In this paper, we used the 
atom-atom interaction Eqs. (|^). This would be modified 
for atoms in a one-dimensional confining potential such 
as an atom waveguide or a cigar-shaped potential. Ac- 
cording to Rcf. , a two-body potential of the atoms 
in such a confining potential can be written as 



2Trh a 1 
m 7r6^ 



1.4603 — 

0| 



(62) 



where b± — ^Jfij muj± . For b± ^ a, which our parameter 
b± ~ 0.29 /im satisfies, Eq. ( |62[ ) is smoothly reduced to 
Eq. d). 



[1] C. J. Myatt, E. A. Burt, R. W. Christ, E. A. Cornell, 
and C. E. Wieman, Phys. Rev. Lett. 78, 586 (1997). 

[2] J. Stenger, S. Inouye, D. M. Stamper-Kurn, H.-J. Mies- 
ner, A. P. Chikkatur, and W. Ketterle, Nature (London) 
396, 345 (1998). 

[3] H.-J. Miesner, D. M. Stamper-Kurn, J. Stenger, S. In- 
ouye, A. P. Chikkatur, and W. Ketterle, Phys. Rev. Lett. 
82, 2228 (1999); D. M. Stamper-Kurn, H.-J. Miesner, A. 
P. Chikkatur, S. Inouye, J. Stenger, and W. Ketterle, 
ibid. 83, 661 (1999). 

[4] T. L. Ho and V. B. Shenoy, Phys. Rev. Lett. 77, 3276 

(1996) . 

[5] B. D. Ersy, C. H. Creene, J. P. Burke, and J. L. Bohn, 

Phys. Rev. Lett. 78, 3594 (1997). 
[6] H. Pu and N. P. Bigelow, Phys. Rev. Lett. 80, 1130 

(1998); 80, 1134 (1998). 
[7] P. Ohberg and S. Stenholm, Phys. Rev. A 57, 1272 

(1998) . 

[8] E. P. Bashkin and A. V. Vagov, Phys. Rev. B 56, 6207 

(1997) . 

[9] E. Timmermans, Phys. Rev. Lett. 81, 5718 (1998). 
[10] B. D. Ersy and C. H. Creene, Phys. Rev. A 59, 1457 

(1999) . 

[11] P. Ao and S. T. Chui, Phys. Rev. A 58, 4836 (1998); 59, 
1473 (1999). 

[12] M. Trippenbach, K. Coral, K. Rzazewski, B. Malomed, 

and Y. B. Band, J. Phys. B 33, 4017 (2000). 
[13] B. Tanatar and K. Erkan, Phys. Rev. A 62, 053601 

(2000) . 

[14] D. Cordon and C. M. Savage, Phys. Rev. A 58, 1440 

(1998) . 

[15] P. Ohberg, Phys. Rev. A 59, 634 (1999). 

[16] C. K. Law, H. Pu, N. P. Bigelow, and J. H. Eberly Phys. 

Rev. Lett. 79, 3105 (1997). 
[17] C. C. Bradley, C. A. Sackett, J. J. ToUett, and R. C. 

Hulet, Phys. Rev. Lett. 75, 1687 (1995). 
[18] H. T. C. Stoof, J. Stat. Phys. 87, 1353 (1997); M. Ueda 

and A. J. Leggett, Phys. Rev. Lett. 80, 1576 (1998); C. 



Huepe, S. Metens, C. Dewel, P. Borckmans, and M. E. 

Brachet, ibid. 82, 1616 (1999); J. A. Freire and D. P. 

Arovas, Phys. Rev. A 59, 1461 (1999). 
[19] M. Olshanii, Phys. Rev. Lett. 81, 938 (1998). 
[20] D. S. Petrov, M. Holzmann, and C. V. Shlyapnikov, 

Phys. Rev. Lett. 84, 2551 (2000). 
[21] E. V. Coldstein, M. C. Moore, H. Pu, and P. Meystre, 

Phys. Rev. Lett. 85, 5030 (2000). 
[22] M. R. Matthews, D. S. Hall, D. S. Jin, J. R. Ensher, C. 

E. Wieman, E. A. Cornell, F. Dalfovo, C. Minniti, and 

S. Stringari, Phys. Rev. Lett. 81, 243 (1998). 
[23] S. Inouye, M. R. Andrews, J. Stenger, H.-J. Meisner, D. 

M. Stamper-Kurn, and W. Ketterle, Nature (London) 

392, 151 (1998); Ph. Courteille, R. S. Freeland, D. J. 

Heinzen, F. A. van Abeelen, and B. J. Verhaar, Phys. 

Rev. Lett. 81, 69 (1998); S. L. Cornish, N. R. Claussen, 

J. L. Roberts, E. A. Cornell, and C. E. Wieman, ibid. 85, 

1795 (2000). 

[24] H. Kleinert, Path Integrals in Quantum Mechanics, 
Statistics, and Polymer Physics, (World Scientific, Sin- 
gapore, 1990). 

[25] In the case of a single condensate, the collective coor- 
dinate space is three dimensional, corresponding to the 
sizes of the condensate along x, y, and z directions. 
One of the authors (Y.Y.) studied the bounce trajec- 
tory by solving the imaginary-time equation of motion 
SSe/SCI = analytically and numerically [^. It was 
shown that the trajectory is approximately straight when 
the number of condensed bosons is slightly below a cer- 
tain critical number. 

[26] F. R. Klinkhamer and N. S. Manton, Phys. Rev. D 30, 
2212 (1984). 

[27] W. Zweger, Z. Phys. B: Condens. Matter 51, 301 (1983). 
[28] Y. Yasui, T. Takaai, and T. Ootsuka, J. Phys. A 34, 2643 
(2001). 

[29] I. Affleck, Phys. Rev. Lett. 46, 388 (1981). 



